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QUASINEUTRAL  PARTICLE  SIMULATION  OF  MAGNETIZED  PLASMA  DISCHARGES: 

GENERAL  FORMALISM  AND  APPLICATION  TO  ECR  DISCHARGES 

1.  Introduction 

In  recent  years,  there  has  been  considerable  interest  in  the  use  of  low  pressure,  high 
density  plasma  sources  for  materials  processing,  e.g.  for  semiconductor  etching,  and  also  for 
some  types  of  chemical  vapor  deposition.*  These  types  of  sources,  which  include  electron 
cyclotron  resonance  (ECR),  helicon  and  inductively-coupled  sources,  operate  at  gas 
pressure  ranging  from  10  mTorr  to  below  1  mTorr,  and  plasma  density  from  several  x  10“ 
to  10*^  cm In  this  regime,  fluid  models  are  not  well  matched  to  the  characteristics  of 
either  the  ions  or  the  electrons.  For  example,  mean  free  paths  can  be  comparable  to  device 
dimensions,  the  electron  response  to  the  driving  electromagnetic  fields  can  be  nonlocal,  and 
velocity  distribution  functions  can  be  non-Max wellian  (particularly  for  the  high  energy 
electrons  that  determine  sheath  potentials,  ionization  rates,  and  other  inelastic  collisional 
processes).  For  these  reasons,  fully  kinetic  simulation  models  are  needed  to  properly 
represent  the  overall  physics  and  chemistry  of  the  discharge. 

A  natural  approach  to  modeling  of  these  types  of  discharges  is  the  use  of  particle-in- 
cell  (PIC)  computer  simulations  which  also  include  a  Monte  Carlo  (MC)  representation  of 
collisional  processes  (PIC/MC  models).^  In  these  simulations,  the  motion  of  particles 
between  collisions  is  followed  deterministically,  under  the  influence  of  specified  external 
magnetic  fields,  and  self-consistent  electrostatic  fields  computed  by  calculating  the  charge 
density  of  the  particles  and  then  solving  Poisson’s  equation.  Recent  years  have  seen  major 
advances  in  the  development  and  use  of  PIC/MC  codes.^'*  *  In  most  cases,  the  objective  of  a 
simulation  is  to  calculate  the  characteristics  of  the  plasma  steady  state,  or  at  most  of  a 
slowly  varying  state,  as  a  function  of  the  machine  control  parameters.  However,  the 
enormous  range  of  spatial  and  temporal  scales  in  a  high-density  discharge  precludes 
straightforward  application  of  the  PIC/MC  procedure  to  simulations  of  the  entire  discharge. 
As  an  example,  the  spatial  and  temporal  scales  for  a  typical  ECR  plasma  are  summarized  in 
Table  1.  With  respect  to  plasma  transport  and  electrical  properties,  the  approach  to  a  steady 
state  occurs  over  time  scales  characterized  by  the  escape  of  ions  to  the  walls,  which  are 
typically  on  the  order  of  0.1  ms  to  several  ms.  The  evolution  to  a  steady  state  of  the 
chemistry  may  occur  on  an  even  longer  time  scale  characterized  by  the  residence  of  neutrals 
Manuscript  approved  June  11,  1997. 
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jjj  the  system.  However,  s  strfti^htforwsrd  PIC  code  must  resolve  &11  of  the  shorter  time 
scales,  down  to  and  including  the  electron  plasma  frequency,  electron  gyrofrequency  and 
microwave  frequency  time  scales,  which  are  on  the  order  of  picoseconds.  Similarly,  the 
spatial  scales  of  real  interest  in  a  reactor-scale  simulation  are  macroscopic  scale  lengths 
(centimeters),  but  a  straightforward  particle  simulation  would  have  to  resolve  a  vast  range  of 
spatial  scales,  down  to  the  electron  gyroradius,  the  Debye  length,  and  the  thickness  of 
passive  sheaths,  which  are  typically  0.01  to  0.1  cm.  It  would  be  essentially  impossible, 
even  with  a  supercomputer,  to  resolve  this  range  of  spatial  and  temporal  scales  within  a 
multi-dimensional  simulation.  To  the  greatest  extent  possible,  one  would  like  to  simulate 
the  macroscopic  features  and  the  approach  to  equilibrium,  while  excluding  the  fast  time 
scales  and  short  spatial  scales  from  the  simulation. 

One  approach  to  the  problem  of  time  scale  diversity  is  the  use  of  implicit  coding  to 
avoid  resolution  of  the  electron  plasma  frequency  time  scale.'^  ’’  However,  PIC  simulation 
techniques  also  face  a  more  fundamental  difficulty  arising  from  the  circumstances  of 
quasineutrality.  The  internal  electrostatic  field  within  the  bulk  plasma  plays  an  essential 
role  in  the  transport  of  plasma  toward  the  walls,  and  therefore  strongly  influences  the 
overall  structure  and  characteristics  of  the  discharge.  This  field  is  typically  less  than  1  V/cm 
in  the  central  plasma,  up  to  a  few  V/cm  in  the  presheath.  When  viewed  in  the  light  of 
Poisson’s  equation. 


=  4ire(ng-nj),  _ 

the  source  of  this  field,  Uj  -  n^,  represents  an  electron-ion  charge  separation  that  is  less  than 
10-5  of  the  density  of  either  electrons  or  ions.  This  type  of  situation  is  actually  quite 
common  throughout  plasma  physics,  as  processes  occurring  on  ion  time  scales  typically  are 
quasineutral,  i.e.  have  I  n^  -  n^  I  «  Uj.  If  the  electrons  and  ions  in  the  discharge  are 
represented  by  simulation  macroparticles,  any  attempt  to  calculate  0  directly  from  Eq.  (1) 
would  be  overwhelmed  by  statistical  noise.  For  example,  in  a  million-particle  2-D 
simulation  with  a  lOQxlOO  grid,  there  are  typically  100  electrons  or  ions  in  each  cell,  and 
therefore  the  statistical  fluctuations  in  nj  -n^  within  a  cell  would  be  on  the  order  of  10%, 
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which  is  four  orders  of  magnitude  larger  than  the  actual  value.  Numerical  schemes 
involving  Poisson’s  equation  are  obviously  very  difficult  (and  actually  inappropriate)  in  the 
quasineutral  limit.  Indeed,  Chen‘*  noted  long  ago  that,  "In  a  plasma,  it  is  usually  possible  to 
assume  n^  =  nj  and  V  •  E  0  at  the  same  time.  This  is  a  fundamental  trait  of  plasmas,  one 
which  is  difficult  for  the  novice  to  understand.  Do  not  use  Poisson’s  equation  to  obtain  E 
unless  it  is  unavoidable!" 

Over  the  years,  this  approach  has  been  followed  in  many  analytic  and  numerical 
models  which  represent  the  plasma  as  a  fluid,  or  represent  the  electrons  as  a  fluid  within 
some  hybrid  scheme.*^  These  methods  circumvent  the  use  of  Poisson  s  equation  by 
neglecting  electron  inertia  and  determining  E  from  the  resulting  simplified  electron 
momentum  conservation  equation.  In  these  models,  nj  is  determined  from  dynamical 
equations,  but  n^.  is  simply  set  equal  to  nj  to  maintain  quasineutrality.  This  procedure 
eliminates  temporal  scales  on  the  order  of  the  electron  plasma  frequency,  as  well  as  spatial 
structures  on  the  Debye  length  scale.  However,  in  kinetic  models,  and  in  particular  in  PIC 
simulations,  the  usual  approach  has  been  to  calculate  ^  from  Poisson  s  equation. 

We  have  developed  a  new  fully  kinetic  approach  to  the  analysis  and  simulation  of 
bounded,  weakly  collisional  plasma  discharges  in  which  the  bulk  plasma  is  quasineutral. 
Our  approach  is  motivated  by  the  quasineutral  fluid  techniques  described  in  the  previous 
paragraph.  At  the  present  stage,  the  model  specifically  treats  axisymmetric  situations  in 
which  the  electrons  are  strongly  magnetized,  as  in  ECR  plasmas,  but  we  believe  the  general 
approach  can  also  be  extended  to  three  dimensions  and  used  for  unmagnetizedat  weakly 
magnetized  cases,  and  is  equally  applicable  to  quasineutral  processes  in  high-temperature 
collisionless  plasma.  Non-uniform  externally  specified  magnetic  fields  are  included,  and 
electrostatic  fields  are  computed  self-consistently,  but  Poisson’s  equation  is  not  used.  In  the 
numerical  implementation,  both  electrons  and  ions  are  represented  as  PIC  particles,  with 
Monte  Carlo  collisions.  The  spatial  gridding  relects  the  macroscopic  scale  lengths,  and  the 
time  steps  are  chosen  to  resolve  particle  motion  over  macroscopic  lengths.  We  shall  also 
give  some  examples  in  which  the  model  provides  a  natural  framework  for  analytic 
calculations,  leading  to  new  insights.  The  model  has  the  following  characteristics. 
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(i)  Electrons  are  treated  as  guiding  center  particles,  and  in  fact  are  firmly  attached  to  a 
single  flux  surface.  This  eliminates  the  electron  gyro  time  scale  from  the  dynamics. 
However,  the  ion  gyro  motion  is  resolved  both  spatially  and  temporally,  and  ions  are  pushed 
through  their  orbits  under  the  influence  of  the  electrostatic  and  Lorentz  forces  (as  well  as 

collisions). 

(ii)  The  electric  field  component  En,  parallel  to  the  magnetic  field  lines,  is  determined 
by  the  requirement  that  E,,  drive  the  electrons  to  maintain  quasineutrality.  This  field  can  be 
specified  by  a  slightly  modified  form  of  the  electron  parallel  momentum  conservation 
equation  with  neglect  of  electron  inertia. 

(iii)  Sheaths,  at  either  grounded  or  floating  surfaces,  are  regarded  as  thin  potential 
barriers  to  electron  flow.  The  sheath  stmcture  is  not  resolved,  but  the  sheath  potentials  are 
determined  self-consistently  by  the  requirement  that  electron  flow  to  walls  be  consistent 
with  preservation  of  quasineutrality  in  the  plasma,  and  also  with  constraints  on  current  flow 

to  the  walls. 

(iv)  The  Bohm  condition  on  ion  flow  at  the  bulk-sheath  interface  is  imposed  as  a 
boundary  condition  on  the  bulk  flow.  This  is  a  key  aspect  of  the  model,  since  the  Bohm 
condition  is  the  principal  constraint  driving  plasma  flow  to  the  walls.  The  quasineutral  pre¬ 
sheath  is  resolved  within  the  model. 

(v) The  relative  plasma  potential  on  different  magnetic  field  lines,  and  therefore  the 
electric  field  component  Ejl  transverse  to  the  magnetic  field,  is  determined  through 
conditions  on  the  transverse  ion  flow  necessary  to  maintain  quasineutrality.  The  transverse 
electric  field  and  the  transverse  ion  flow  depend  significantly  on  whether  the  field  lines 
terminate  on  insulating  or  conducting  walls.  In  the  case  of  conducting  walls,  the  sheath 
potentials  and  the  transverse  electric  field  are  intimately  related  and  must  be  determined 
self-consistently.  In  this  paper,  we  present  the  formalism  for  both  cases,  and  show 
simulation  results  for  the  insulating  case,  which  is  considerably  simpler  to  implement. 

Resonant  heating  of  the  plasma  by  the  microwaves  can  also  be  represented  in  a 
linearized  formulation  which  integrates  analytically  over  the  microwave  and  electron  gyro 
frequencies.  However,  this  aspect  is  not  discussed  in  the  present  paper. 
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Within  this  framework,  plasma  oscillations  are  excluded  from  the  model.  Simulation 
time  steps  need  not  resolve  the  electron  plasma  or  gyro  frequency;  thus  they  can  be  chosen 
to  resolve  motion  of  individual  particles  over  macroscopic  lengths  of  interest,  and/or 
collisional  time  scales  of  interest.  In  typical  cases,  we  use  ion  time  steps  of  order  few  x  10  ’ 
sec,  and  sub-cycled  electron  time  steps  of  order  lO  ®  sec.  The  result  is  a  very  efficient  code 
that  can  be  run  to  times  of  real  interest,  e.g.  milliseconds,  in  a  few  hours  on  a  workstation. 
The  self-consistent  separation  of  sheath  and  bulk  plasma  also  facilitates  intuitive 
understanding  and  provides  a  useful  framework  for  analytic  calculations.  For  example,  m 
Sec.  7  we  shall  see  how  self-consistent  variations  of  the  sheath  potential  can  strongly  inhibit 

cross-field  ion  transport. 


2,  Structure  and  Geometry  of  the  Model 

To  help  motivate  the  assumptions  of  the  model,  it  may  be  useful  to  consider  some 
typical  plasma  conditions  which  are  of  interest  to  us  (while  keeping  in  mind  that  the  model 
is  applicable  to  a  much  broader  range  of  conditions).  In  ECR  processing  sources,  the 
microwave  frequency  is  usually  2.45  GHz,  so  that  electron  gyroresonance  occurs  at  875  G. 
In  a  two-solenoid  source  (Fig.  1),  the  magnetic  field  lines  typically  diverge  downstream,  so 
that  1000  G  where  the  microwaves  are  introduced,  decreases  to  the  resonant  value 
nearby,  and  falls  to  as  low  as  20  G  far  downstream.  The  electron  temperature  T^  is  typically 
a  few  eV.  The  ion  temperature  T^  is  typically  much  lower,  but  the  ion  flow  speed  reaches 
the  ion  sound  speed  c,  =  (T^mi)^'’  at  the  interface  between  the  bulk  plasma  and  the  sheath. 
Many  gas  compositions  are  used  for  processing,  and  pressure  may  vary  from  0.1  mTorr  to  > 
10  mTorr,  but  for  specificity  we  shall  consider  Ar  at  pressure  1  mTorr.  Typically,  the 
plasma  density  will  be  few  x  10"  to  \ so  that  the  ionization  fraction  is  between  1% 

and  20%. 

Under  these  conditions,  the  electron  gyroradius  is  extremely  small  (<  1  mm)  and  the 
electron  mean  free  path  is  so  long  (about  40  cm  for  electron-neutral  collisions  and  30  cm 
for  electron-ion  collisions  at  plasma  density  iC’cm'^  and  T,=4eV)  that  electron  collisional 
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diffusion  across  field  lines  is  very  slow.  Typically  an  electron  will  diffuse  across  field  lines 
by  <  1  cm,  before  eventually  escaping  from  the  system  to  a  wall.  Furthermore,  we  shall 
restrict  our  attention  in  this  paper  to  axisymmetric  systems,  so  that  all  collisionless  electron 
drifts  are  azimuthal.  Thus,  it  is  quite  reasonable  to  assume  that  each  electron  is  strictly 
confined  to  a  particular  flux  tube.^^  In  the  framework  of  r-z  geometry,  where  the  azimuthal 
coordinate  is  ignorable,  we  may  think  of  the  electrons  as  strung  out  along  the  magnetic  field 
lines.  Thus  the  electron  motion  is  essentially  one-dimensional  in  the  r-z  plane. 

In  order  to  treat  the  electron  dynamics  efficiently  and  accurately,  we  use  magnetic 
field  lines,  denoted  by  a  discrete  index],  as  one  of  the  coordinates  for  our  grid.  The  spacing 
of  the  field  lines  chosen  for  the  grid  is  arbitrary,  and  in  principle  may  be  chosen  to  optimize 
the  resolution  in  regions  of  particular  interest.  However,  there  are  also  numerical 
constraints  on  the  choice  of  grid  field  lines,  e.g.,  it  is  necessary  to  maintain  an  adequate 
number  of  particles  on  each  grid  field  line  to  control  statistical  fluctuations.  The  grid  used 
in  typical  simulations  is  shown  in  Fig.  2.  The  axial  coordinate  z  is  used  as  the  other 
coordinate  for  the  grid;  equally  spaced  grid  lines  in  z  are  denoted  by  the  index  k.  The  use 
of  a  non-uniform  and  non-orthogonal  grid,  which  also  intercepts  the  radial  wall  at  oblique 
angles,  complicates  the  code  structure  considerably,  but  affords  great  simplicity  and 
efficiency  in  the  treatment  of  electron  dynamics:  each  electron  is  permanently  associated 
with  a  single  field  line  j  during  its  entire  lifetime  in  the  system,  and  in  the  particle  lay-down 
is  linearly  distributed  onto  the  two  nearest  grid  points,  (j.k)  and  (j,k+l). 

The  ion  gyroradius  ranges  from  a  few  mm  for  thermal  ions  in  the  strong  field  region, 
to  tens  of  cm  for  ions  moving  at  c*  in  the  weak  field  region.  Furthermore,  the  ion  mean  free 
path  for  charge  exchange^^’^'^  (  ~  10  cm)  is  often  comparable  to  or  moderately  larger  than  the 
ion  gyroradius.  Thus  the  ions  are  only  weakly  magnetized.  In  the  model,  the  ion  orbits  are 
computed  in  full  3-D,  under  the  influence  of  the  externally  specified  magnetic  field,  the  self- 
consistent  electric  fields  which  are  interpolated  from  the  field  line  grid  onto  an  r-z  grid,  and 
Monte  Carlo  collisions.  After  the  completion  of  an  ion  time  step,  the  ion  density  is  laid 
down  on  the  j-k  grid,  as  follows.  Each  ion  is  attributed  in  z  fully  to  its  nearest  neighbor  grid 
value  k,  and  then  is  apportioned  among  the  two  field  lines  between  which  it  is  located, 
according  to  the  quadratic  formula 
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(2a) 


-  r: 


jk 


_  r^  -  rjk  (2b) 

^j+l,k  “2  _  j~2  ' 

^j+l,k  ^jk 

which  is  appropriate  to  r-z  geometry.  Here,  the  ion  coordinate  r  is  located  between  fieid 
iines  j  and  j+1,  with  r,,  and  tj„*  the  coordinates  of  the  nearest  grid  points  on  fieid  iines  j  and 
j+1  respectiveiy,  and  Pji  is  the  fracUon  of  the  ion  that  is  attributed  to  field  line  j. 


3.  Electron  Dynamics  and  Parallel  Electric  Field 

Since  the  electrons  are  represented  as  guiding  center  particles,  each  electron  is 
characterized  by  its  location  along  the  field  line.  In  the  code  we  use  the  Cartesian 
coordinate  z  to  specify  this  location,  but  in  analytic  development  it  is  often  useful  to  use  $, 
the  curvilinear  coordinate  along  the  field  line.  In  addition,  each  electron  has  a  parallel 
velocity  v„  and  a  perpendicular  (gyrating)  velocity  v^.  We  assume  that  the  magnetic 
moment  m  =  12 1 B I  is  a  constant  of  the  motion  in  between  collisions,  so  usually  it  is 

convenient  to  characterize  the  electron  by  its  value  of  m  rather  than  Vj^.  The  effective 
parallel  force  acting  on  an  electron,  between  collisions,  is  then 

F|,  =  -eE„-M(dB/d$), 

where  the  second  term  is  the  mirror  force.  In  the  computational  model,  the  electrons  are 
pushed,  between  collisions,  as  simulation  particles  subject  to  the  force  F,,.  The  time  step  At* 
for  the  electron  push  is  chosen  to  resolve  electron  motion  over  macroscopic  lengths  of 
interest,  and  also  to  resolve  electron  collisional  time  scales.  Typically  At^  is  a  small  fraction 
of  the  time  step  Atj  used  to  push  ions. 

Let  us  assume  for  the  moment  that  the  number  of  electrons  on  a  field  line  j  is  equal  to 
the  number  of  ions  assigned  to  that  field  line,  so  that  the  field  line  is  globally  quasineutral. 
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(We  will  later  discuss  how  to  insure  that  this  is  true.)  The  electron  Debye  length  is  always 
smaller  than  any  scale  length  resolved  in  the  model,  and  the  electron  plasma  frequency  is 
always  fast  compared  to  any  time  scale  resolved  in  the  model.  Thus,  if  there  were  any 
departure  from  local  quasineutrality,  the  resulting  strong  electric  field  would  drive  electron 
currents  parallel  to  B,  which  would  restore  quasineutrality  within  a  time  scale  of  several 
electron  plasma  periods,  i.e.  essentially  instantaneously  on  the  time  scale  of  the  model. 
Thus,  the  macroscopic  parallel  electric  field  always  takes  the  value  necessary  to  keep  the 
electron  density  n^  equal  to  the  ion  density  n^  To  specify  this  electric  field,  we  can  begin 
with  the  electron  momentum  equation.  For  magnetized  electrons  on  curved  field  lines,  this 
equation  takes  the  form 


-eEii  = 


IB  I  a  (P^H+ngmgO^  - 


1  a 


+ -^  (riemeUeii)  -  (4) 


Hg  a  s  I B I 

Here,  m^  is  the  electron  mass,  PgH  =  ngTeH  is  the  electron  parallel  pressure. 


Pelt  =  X dV|,me(V|rUell)^  fe(V||.Vx), 


(5) 


Ugii  is  the  mean  electron  fluid  velocity,  m  is  the  average  magnetic  moment  at  the  specified 
location,  and  is  the  mean  electron  momentum  transfer  collision  frequency. 

Within  the  particle  simulation,  the  first  three  terms  on  the  right  hand  side  (RHS)  of 
Eq.  (4)  can  be  evaluated  at  each  point  of  the  grid,  by  laying  down  the  mean  quantities  for 
the  electrons  assigned  to  that  grid  point.  These  terms  represent  the  "ambipolar"  electric 
field  -eEi*!®^  necessary  to  sustain  the  quasineutrality  relation 


if  Uj  were  a  specified  time-independent  function  of  $ .  The  remaining  term  (the  inertial  term) 
is  a  small  correction,  of  order  mg/m;,  where  mj  is  the  ion  mass.  In  quasineutral  fluid 
formulations,  where  n^  is  simply  set  equal  to  Uj,  the  inertial  term  is  neglected.  However  this 
is  unsatisfactory  in  a  particle  simulation,  where  n^  and  nj  are  separately  determined  by  the 
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particle  evolution;  if  the  inertial  term  is  neglected,  n,.  and  nj  will  drift  apart  as  nj  slowly 
changes  and  n^  does  not  respond.  Even  worse,  in  a  particle  code  both  n,.  and  tij  are  subject  to 
continual  fluctuations,  and  the  ambipolar  field  allows  n^  and  nj  to  separate  through  these 
fluctuations. 

There  are  several  possible  approaches  to  the  calculation  of  the  inertial  term,  or  of 
some  approximate  form  that  couples  n^  to  n^.  In  a  previous  publication,^^  we  developed  an 
approximate  method  in  which  the  electrons  are  pushed  in  the  ambipolar  field,  and  then  at 
the  end  of  the  time  step  a  correction  field  is  applied  that  is  chosen  to  restore  quasineutrality. 
Although  this  method  worked  well  in  most  respects,  we  have  found  another  approach  that  is 
even  simpler,  and  very  accurately  conserves  energy  over  long  time  scales.  We  simply  drop 
the  inertial  term  in  (4)  and  substitute  tij  for  n^  in  the  first  term  of  (4),  so  that 


IB  I  9  [n,  (T„ii+m^Ueii^]  .  7.  9  I B I 
-eE  II  =  - •=-=•  - M  ST- 
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IBI 
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VenieUeii . 


(7) 


Equation  (7)  keeps  the  electron  density  closely  coupled  to  the  ion  density,  while  avoiding 
the  very  high-frequency  plasma  oscillations  that  are  linked  to  the  last  term  of  (4).  This  can 
be  seen  by  substituting  (7)  for  En  in  the  exact  momentum  equation  (4),  to  obtain 


in„ 


at 


n, 


(8) 


Equation  (8)  shows  that  the  electrons  are  always  accelerated  up  the  gradient  in  n/n^,  i.e. 
toward  the  point  of  maximum  positive  net  charge  density,  thereby  restoring  quasineutrality. 
We  showed  in  Ref.  26  that  this  scheme  is  stable,  preserves  while  supporting  only 
lower-frequency  oscUlations,  conserves  energy  to  good  accuracy,  and  accurately  represents 
kinetic  effects.  For  example,  the  usual  dispersion  relation  for  ion  sound  waves,  including 
the  Landau  damping  terms,  can  be  derived.^  We  also  gave  several  examples  of  successful 


implementation  of  the  scheme  in  a  one-dimensional  particle  simulation  with  no  boundaries, 

26 

including  linear  and  nonlinear  ion  sound  waves,  and  free  expansion  of  a  plasma. 
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4.  Sheaths  at  Passive  Surfaces 


At  any  material  surface  in  contact  with  the  plasma,  there  is  a  sheath  where 
quasineutrality  fails.  We  shall  consider  here  only  sheaths  at  unpowered  or  weakly  dc-biased 
surfaces,  e.g.  vessel  walls.  The  sheath  thickness  is  typically  a  few  times  the  Debye  length 
Xd,  which  is  very  thin  for  conditions  of  interest.  For  example,  Xp  =  15  Mtn  at  =  4  eV, 
n  =  10‘2  cm■^  Such  a  sheath  can  be  treated  as  simply  a  thin  potential  jump  <t>^,  which 
accelerates  positive  ions  toward  the  wall,  but  reflects  most  electrons  back  into  the  plasma. 
At  any  given  time,  only  a  few  electrons  with  high  enough  parallel  energy  can  surmount  the 
sheath  potential  barrier  and  reach  the  wall.  (However,  any  electron  can  eventually  diffuse 
up  to  high  energy,  due  to  ECR  heating  or  electron-electron  collisions,  and  escape  from  the 
system.)  The  effect  of  the  sheath  on  the  plasma,  or  on  the  surface,  is  essentially  completely 
characterized  by  the  sheath  potential. 

To  determine  the  value  of  the  sheath  potential,  we  use  an  elaboration  of  the  "logical 
sheath"  scheme  of  Parker,  Procassini  and  Birdsall.^^  Our  formulation  incorporates  the 
Bohm  flow  criterion,  applies  appropriately  to  both  conducting  and  insulating  walls,  and  can 
be  used  for  a  multi-dimensional  plasma  with  magnetized  electrons.  The  basic  idea  is  that 
the  sheath  potential  takes  the  value  which  allows  the  "correct"  electron  flux  to  reach  the 
wall.  The  conditions  for  determining  the  correct  flux  depend  on  whether  the  wall  is 
conducting  or  insulating. 

A.  Insulating  walls 

An  insulating  wall  exposed  to  a  plasma  acquires  a  surface  charge  during  an  initial 
transient  period  (which  we  do  not  resolve),  and  thereafter  the  electrical  current  density  to 
any  point  on  the  wall  must  be  zero.  Thus  the  sheath  potential  at  each  end  of  field  line  j 
must  take  the  value  (the  floating  potential)  which  sets  the  flux  of  electrons,  through  the 
sheath  to  the  wall,  equal  to  the  ion  flux  into  the  sheath  from  the  plasma.  (All  positive  ions 
which  reach  the  boundary  of  the  simulation  pass  through  the  sheath  and  reach  the  wall.) 
Furthermore,  according  to  the  Bohm  criterion,^*  ions  must  flow  from  the  bulk  plasma  into 
the  sheath  at  a  mean  velocity  equal  to  the  ion  sound  speed  c^,  i.e.  the  ion  flux  to  any  wall  is 
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j  Bohm  _ 


(9) 


In  the  next  section  we  shall  discuss  the  way  in  which  the  Bohm  criterion  is  imposed  as  a 
boundary  condition  on  the  ion  flow  in  the  model.  However,  the  implication  for  the 
electrons  is  that  the  electron  flux  to  any  wall,  must  satisfy 

00 

=  sin  e  I  dV|(V||fe(V||)  = 

where 


V,  =  (2e0/irg‘^, 

and  0  is  the  angle  between  the  magnetic  field  and  the  wall.  It  the  electron  velocity 
distribution  funcUon  4(v„)  at  the  end  of  field  line  j  is  known,  Eqs.  (10)  and  (1 1)  determine 
the  sheath  potential 

To  implement  this  prescription  numerically,  we  first  push  both  the  electrons  and  ions 
through  a  complete  ion  time  step  Aq.  If  an  electron  reaches  a  wall,  it  is  reflected  back  along 
its  field  line,  on  the  tentative  assumption  that  its  energy  was  insufficient  to  penetrate  the 
sheath  potential  barrier.  However,  it  is  noted  that  that  particular  electron  belongs  to  the  set 
S  bounce  of  electrons  on  field  line  j  that  reflected  off  the  wall  during  Atj.  In  general,  this 
number  of  electrons  will  be  much  larger  than  the  number  of  ions  that  reach  that  wall  during 
the  time  step,  AN^^^""  h  Ji®°'“"AtiAj,  where  Aj  is  the  area  intercepted  on  the  wall  by  the  flux 
tube  associated  with  field  line  j.  We  then  identify  the  subset  of  consisting  of 

the  ANj®*'"’  electrons  which  had  the  highest  parallel  energy  at  the  moment  of  reflection. 
These  electrons  are  assumed  to  have  escaped  during  the  time  step;  hence  they  are  discarded. 
The  sheath  potential  is  set  equal  to  the  lowest  parallel  kinetic  energy  of  any  of  the 
electrons  in  the  subset  ,  at  the  time  it  reached  the  simulation  boundary. 

There  are  a  number  of  numerical  points  which  must  be  considered  in  actually 
implementing  this  scheme.  The  sheath  potential  depends  on  the  high-energy  tail  of  the 
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electron  distribution,  which  always  contains  relatively  few  simulation  electrons  and  is 
susceptible  to  noise.  To  minimize  this,  Atj  must  be  large  enough  that  a  substantial  absolute 
number  of  electrons  are  allowed  to  escape  during  the  time  step.  However,  the  fraction  of 
electrons  which  leave  during  At;  must  be  small.  In  a  steady  state  situation,  the  departed 
electrons  will  be  replaced  by  new  electrons  generated  by  ionizing  collisions,  and  one  must 
be  careful  not  to  overcount  ionizing  collisions  involving  electrons  that  actually  escaped 
during  Atj.  Also,  there  may  be  electrons  with  high  parallel  velocities  that  bounce  off  both 
walls  during  the  time  step  Atj,  and  these  must  be  accounted  for  correctly. 


B.  Conducting  walls 

When  the  walls  are  conducting,  there  will  typically  be  a  boundary  condition 
* 

specifying  the  potential  on  any  given  segment  of  the  wall.  In  the  simplest  case,  the  walls 
may  all  be  regarded  as  grounded,  but  in  other  cases  a  dc  bias  may  be  applied.  (We  do  not 
consider  ac  biases  in  this  paper.)  For  conducting  walls,  there  is  a  more  complicated  set  of 
conditions  that  determines  the  sheath  potentials.  It  is  possible  for  non-zero  electrical  current 
density  to  flow  to  any  given  point  on  the  wall,  i.e.  the  electron  flux  to  the  wall  need  not  be 
locally  equal  to  the  ion  flux  out  of  the  plasma.  Within  the  plasma,  there  may  be  currents 
flowing  along  a  field  line,  and  also  ion  currents  flowing  across  field  lines.  However, 
quasineutrality  imposes  a  constraint  on  the  electron  and  ion  fluxes:  the  total  number  of  ions 
and  electrons  on  any  given  field  line  must  remain  equal.  This  is  the  condition  that 
determines  the  sheath  potentials. 

Let  Jej  “‘(0)  be  the  electron  flux  to  the  wall  at  the  z = 0  end  of  field  line  j,  and  Jej"‘(Zmax ) 
be  the  flux  to  the  wall  at  the  other  end  of  field  line  j.  Let  and  Jij“‘“"(Zn,ax)  the 

Bohm  flux  of  ions  into  the  sheaths  at  z  =  0  and  z^^^  respectively.  Then  maintenance  of 
global  quasineutrality  on  field  line  j  requires  that 


J°f  (0)  +  Jif  = 


/■^max 


Jf  9*^(0)  +  (Zn,ax)  + 


dz  V- Jixj  (z)  . 
0 


(12) 
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The  last  term  in  Eq.  (12)  represents  the  net  decrease  in  the  number  of  ions  attributed  to  field 
line  j,  due  to  ion  flow  across  field  lines  within  the  plasma. 

In  Eq.  (12),  the  electron  fluxes  to  the  walls  are  determined  by  the  sheath  potentials 

^sj(0)  and  ^sj(Zmax)* 


.-Vsj(z=0) 


j°Vt(z=0)  =  sin 


jo 


dV||V||fe(z=0,  V||)  , 


(13a) 


00 


Jef  (Zmax)  =  sin 


dV||V,|fe(Z^,V||)  , 

(  ^max) 


(13b) 


where 


Vsj(z=0)  =  [2c4>y(z=0)/ 

Vsj(Zmax)  = 


0JQ  is  the  angle  between  field  line  j  and  the  wall  at  z  —  0,  and  max  angle  at  z  -  z^ax-  If 
the  walls  are  grounded  («=0),  then  the  plasma  potential  just  inside  the  sheath,  at  either  z  =  0 
or  Zmax.  is  just  equal  to  the  sheath  potential.  However,  the  difference  between  the  plasma 
potentials  0j(Zmax)and  ^j(O)  is  given  by  the  line  integral  of  En(z),  the  parallel  electric  field 
given  by  Eq.  (7).  Thus  the  two  sheath  potentials  are  related  by 


^sj(Zmax)-^sj(0) 


•^max 

dzE  II  (z)  , 

•'o 


(14) 


and  Eq.  (12)  in  effect  specifies  both  sheath  potentials,  provided  the  electron  energy 

distributions  and  the  ion  flux  are  known  from  the  simulation. 

If  there  is  an  externally  imposed  potential  difference  along  field  line  j,  between  the 
walls  at  z  =  0  and  z^ax.  then  Eq.  (14)  takes  the  more  general  form 
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=  ^ext 


■^max 

dzE||  (z)  , 

■*0 


(15) 


where  is  the  imposed  potential  difference.  While  Eq.  (12)  specifies  the  total  electron 
flux  escaping  from  the  field  line,  Eq.  (15),  together  with  the  rest  of  the  simulation  equations, 
determines  separately  the  flux  of  electrons  to  the  walls  at  z  =  0  and  and  thus  in  effect 
the  electric  current  flowing  along  the  field  line.  In  this  way,  the  simulation  calculates  the 
electric  current  density  in  terms  of  the  imposed  voltage,  i.e.  it  gives  a  macroscopic  Ohm’s 
law  for  the  plasma. 

It  is  clear  from  this  formulation  that  in  the  case  of  conducting  walls  the  sheath 
potentials  are  intimately  connected  to  the  cross-field  ion  flow.  We  shall  see  that  it  is 
necessary  to  calculate  both  the  sheath  potentials  and  the  ion  flow  self-consistently  at  each 
time  step.  Therefore  we  shall  defer  the  discussion  of  a  scheme  for  actually  calculating  the 
sheath  potentials  to  Sec.  6,  where  the  ion  dynamics  are  discussed. 


5.  Imposition  of  the  Bohm  Criterion  as  a  Boundary  Condition 

In  Sec.  4,  we  have  shown  how  to  specify  the  correct  sheath  potential,  and  impose  this 
potential  barrier  as  a  boundary  condition  to  the  electron  flow.  In  a  real  plasma,  or  in  a 
model  that  calculates  the  internal  electric  field  through  Poisson’s  equation,  a  sheath  with  the 
correct  potential  will  very  quickly  form,  as  a  response  to  the  rapid  initial  loss  of  electrons  to 
the  wall.  Simultaneously,  the  ions  will  be  accelerated  to  a  mean  speed  u>,  normal  to  the 
wall,  given  by 


at  the  interface  between  the  quasineutral  bulk  plasma  and  the  ion-rich  sheath.  Indeed,  it  is 
well  known  that  a  steady-state  monotonic  sheath,  with  n;  >  n^,  can  form  only  where  the  ion 
flow  in  the  bulk  plasma  satisfies  Eq.  (16).  This  is  known  as  the  Bohm  criterion. 
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In  a  magnetized  plasma,  an  ion-rich  sheath  will  form  at  every  surface  which  the 
magnetic  field  lines  intersect  at  an  angle  greater  than  (m^/  As  long  as  there  is  an  ion- 
rich  sheath,  the  Bohm  criterion  (16)  normally  holds,  independently  of  the  direction  of  the 
magnetic  field;  in  (16),  we  use  the  notation  uA  with  superscript  J.  to  denote  the  velocity 

normal  to  the  wall,  not  to  the  magnetic  field. 

However,  our  quasineutral  model,  as  described  in  the  preceding  sections,  describes 
only  the  bulk  plasma  and  simply  assumes  the  existence  of  a  sheath,  without  requiring  that 
the  appropriate  bulk  flow  conditions  obtain.  Equation  (16)  is  not  a  consequence  of  the 
quasineutral  model;  steady  state  solutions  exist  which  satisfy  Eq.  (16),  but  there  are  also 
solutions  which  exhibit  other  ion  flow  patterns  at  the  walls.  For  example,  if  the  ions  are 
initially  cold  and  uniform  and  the  electrons  are  isothermal  throughout  the  bulk  plasma,  the 
electron  pressure  will  remain  uniform,  no  electric  field  will  be  generated,  and  there  will  be 
no  ion  flow  toward  the  wall.  Thus,  within  our  model,  the  Bohm  condition  must  be  imposed 
externally  as  a  boundary  condition  to  the  ion  flow,  which  is  necessary  to  specify  the  correct 
solution  of  the  equations.  One  may  say,  from  the  perspective  of  the  bulk  plasma,  that  the 
Bohm  condition  causes  the  flow  of  the  plasma  toward  the  walls. 

Within  the  context  of  a  fluid  model  of  the  quasineutral  bulk  plasma,  Eq.  (16)  can  be 
imposed  as  a  boundary  condition  to  the  ion  momentum  equation.  But  how  is  one  to  impose 
Eq.  (16)  as  a  boundary  condition  within  a  particle  simulation?  We  have  used  a  procedure 
that  seems  to  imitate  the  dynamical  process  that  sustains  the  Bohm  condition  within  the 
quasineutral  bulk  part  of  a  real  discharge. 

In  a  real  plasma,  the  ion  flow  speed  is  equal  to  Cj.  at  the  edge  of  the  bulk  plasma,  but 
falls  off  to  a  much  lower  value  within  the  interior  of  the  plasma.  The  continuity  equation 
then  implies  that  n^  falls  off  at  the  edge  of  the  bulk  plasma.  Quasineutrality  requires  that  n^ 
=  Uj  =  n.  Assuming  T^  does  not  vary  as  rapidly  as  n  at  the  plasma  edge,  the  fall-off  in 
electron  pressure  at  the  edge  leads  to  an  ambipolar  electric  field  [first  term  of  Eq.  (7)]  which 
sustains  the  ion  flow  in  the  "presheath"  region.  In  the  model,  we  define  a  thin  strip  of  depth 
d,  adjacent  to  any  material  surface,  with  d  chosen  to  be  macroscopically  small  but  larger 
than  CjAtj.  Within  the  strip,  we  lay  down  the  ion  velocities  to  determine  the  average  ion 
flow  speed  normal  to  the  surface,  on  field  line  j.  We  then  add  an  increment  8Uif  to  all  of 
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the  ion  velocities,  such  as  to  bring  Uij-  to  c^.  This  procedure  insures  that  Uj  falls  off  at  the 
perimeter  of  the  bulk  plasma  just  as  it  would  if  Bohm  flow  were  occurring.  The  resulting 
pressure  drop  at  the  perimeter  leads  (over  a  time  scale  characterized  by  the  ratio  of  the 
presheath  width  to  c^)  to  the  formation  of  an  appropriate  presheath.  Bohm  flow  is  thereafter 
self-sustaining,  but  only  neutrally  stable,  with  the  increment  Su-J-  a  small  adjustment  at  each 
time  step. 

The  Bohm  flow  condition,  together  with  the  sheath  potential  as  specified  in  Sec.  4, 
together  represent  a  complete  set  of  boundary  conditions  for  the  quasineutral  bulk  plasma. 
We  shall  refer  to  this  as  the  "Bohm  logical  sheath"  procedure. 


6.  Transverse  Electric  Field,  and  Ion  Transport  Across  Magnetic  Field  Lines 

In  Sec.  3  we  showed  how  to  determine  the  parallel  electric  field  En  from  the 
requirements  of  quasineutrality,  as  applied  to  the  electron  flow  along  field  lines.  Thus,  the 
relative  potential  between  any  two  points  on  the  same  field  line  is  known.  If  the  absolute 
potential  were  specified  at  any  single  reference  point  on  each  field  line,  say  at  z  =  0,  then  it 
would  be  known  everywhere,  and  the  transverse  electric  field  Ej^  could  also  be  computed. 
This  additional  information  can  also  be  obtained  from  the  requirements  of  quasineutrality, 
but  it  is  bound  up  with  the  ion  flow,  rather  than  the  electron  flow .  The  conclusions  also 
depend  fundamentally  on  whether  the  field  lines  end  on  insulating  walls  or  on  grounded 
conducting  walls. 

A.  Insulating  walls 

In  the  case  of  insulating  walls,  electrical  current  cannot  flow  from  the  plasma  to  any 
point  on  the  walls.  Furthermore,  in  our  model  no  transverse  electron  current  is  permitted 
within  the  plasma.  It  follows  that  if  there  is  non-zero  net  ion  flux  onto  any  field  line  within 
the  plasma,  a  violation  of  quasineutrality  will  occur  on  that  field  line.  In  a  real  plasma,  a 
large  potential  would  immediately  build  up  all  along  the  field  line,  leading  to  a  reversal  of 
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the  net  transverse  ion  flow  onto  the  field  line,  and  subsequently  to  a  train  of  ion  plasma 
oscillations.  To  model  this  as  an  effect  within  a  quasineutral  framework,  we  shall  stipulate 
that  the  integrated  net  ion  flux  onto  any  field  line  is  zero,  i.e. 


/dzVx- (nUi)  =  0.  '' 

Equation  (17)  is  well  satisfied  on  time  scales  longer  than  the  ion  plasma  frequency, 
expresses  the  condition  of  quasineutrality,  and  can  be  used  to  solve  for  the  reference 
potential,  ^(z=0),  on  each  field  line. 

To  implement  Eq.  (17)  numerically,  we  begin  with  the  known  electric  field  at  the 
previous  time  step  t"*,  assume  that  quasineutrality  was  satisfied  at  t"*,  and  that  the  flux 
condition  (17)  was  satisfied  at  the  previous  half-time  step  First,  we  push  each  ion  to  a 
new  velocity  v*  at  time  using  the  electric  fields  determined  at  t*",  and  to  a  new 

position  (r*,z*)  at  time  At  the  end  of  the  time  step,  there  may  be  a  small  inequality  of 
order  Atj^  between  the  number  of  ions  Ny*  and  electrons  N^j  on  field  line  j,  since  the  flux 
continuity  condition  is  satisfied  exactly  at  time  but  may  be  inaccurate  by  order  Atj 
during  the  time' step: 

ANj*  =  Ny*  -  Nej*  0  (order  At;^). 

This  is  an  indication  that  a  correction  potential  6^j  should  have  been  added  to  the  assumed 
potential  on  field  line  j,  so  as  to  maintain  the  flux  continuity  equation  (17)  throughout  the 
time  step.  This  may  be  thought  of  as  the  correction  to  the  reference  potential  <ij(z=0).  The 
corresponding  correction  to  the  radial  electric  field  is 

6E,jk  =  (5^j-  «^j+i)/(rj+i.k-rjk) 

at  z-grid  point  k,  between  field  lines  j  and  j+1.  This  in  turn  leads  to  a  correction  velocity 


5\jk  =  (e5E,jk/mi)Ati 


(20) 
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and  a  position  displacement 


6r^,  =  (e6E,j,/mi)Ati2  (21) 

for  every  ion  in  this  grid  cell.  Using  Eqs.  (19)  -  (21),  the  resulting  correction  to  the  number 
of  ions  on  field  line  j  can  then  be  written  as 


6Nijj  =  -  (aj+aj.,)6^j  +  aj60j+j. 


(22) 


where  the  coefficients  aj  are  determined  by  first  summing  up  the  displacements  of  the 
particles,  enumerated  with  the  index  n,  which  are  attributed  to  z-grid  point  k  and  lie  between 
field  lines  j  and  (j+1),  and  then  summing  up  the  contributions  of  each  grid  point  k  to  the 
total  number  of  ions  on  field  line  j.  Given  that  the  ions  are  apportioned  among  adjacent 
field  lines  according  to  the  quadratic  scheme  (2),  the  coefficients  are  found  to  be^^ 


The  quantity  SNjjj 
condition  becomes 


__ - 2 - - - i  ^ 

from  Eq.  (22)  is  then  set  equal  to  -ANj*,  so  that  the  quasineutrality 


2r„+ 


dz 


■HSU- 


(23! 


-(aj  +  aj.i)6^j  +  aj8^i --ANj .  (  ) 

Equation  (24)  is  a  linear  ordinary  difference  equation  which  determines  S0j,  the  correction 
to  the  potential  ^j(z  =  0)  which  is  necessary  to  maintain  global  quasineutrality In  this 
way,  the  potential  difference  between  field  lines  is  specified.  After  60j  has  been  calculated, 
the  ions  are  pushed  again  in  the  corrected  fields.  Each  field  line  j  will  then  have  Ny  —  N^j  at 
time  and  Eq.  (17)  will  be  satisfied  at 

It  will  be  noted  that  this  procedure  determines  the  plasma  potentials  and  the  ion  flow 
in  a  coupled  self-consistent  way.  However,  the  sheath  potentials  and  the  electron  flow  play 
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no  role  in  this  deternaination.  As  we  shall  sec,  the  situation  is  different  when  the  walls  of 
the  vessel  are  conducting;  in  this  case,  the  sheath  potentials  and  the  escape  of  electrons  to 
the  walls  must  also  be  included  self-consistently  in  the  calculation. 

B.  Grounded  walls 

In  the  case  of  grounded  walls,  the  wall  serves  as  a  zero  point  for  the  potential  that  is 
common  to  all  field  lines.  The  potential  in  the  bulk  plasma  on  field  line  j,  just  inside  the 
sheath  at  z=0,  is  thus  equal  to  the  sheath  potential  If  is  known,  the  potential  can  be 
calculated  at  any  point  on  the  field  line,  and  the  transverse  E-field  can  then  be  calculated  by 
taking  differences  of  potentials  on  adjacent  field  lines.  However,  the  sheath  potential  is 
determined  from  Eqs.  (12)  -  (14),  which  involve  the  ion  transverse  flow.  Thus  the 
transverse  ion  flow  is  coupled  to  the  escape  of  electrons  through  the  sheaths. 

To  numerically  implement  the  solution  to  this  coupled  problem,  let  us  begin  by 
rewriting  Eq.  (12)  in  the  discrete  form 

AN4"‘  =  ANf  -  ANjjy, 

Here  AN  ®“*  is  the  total  number  of  electrons  on  field  line  j  that  escape  to  the  walls,  through 
the  sheaths  at  both  ends,  during  time  Aq.  For  any  given  electron  energy  distribution,  AN^f  * 
is  a  strong  nonlinear  function  of  ANf is  the  total  number^  ions  that  escape  to  the 
walls  from  flux  tube  j  during  time  Atj.  ANjj^  is  the  net  change  in  the  number  of  ions 
attributed  to  field  line  j,  due  to  transverse  ion  flow  within  the  plasma. 

The  numerical  solution  of  Eq.  (25)  is  similar  in  spirit  to  the  approach  used  for  the  case 
of  insulating  walls,  but  much  more  complicated  in  practice.  Let  us  assume  that  the  plasma 
was  quasineutral  at  the  previous  time  step  t™,  that  the  sheath  potentials  (l>^  are  known  at  t  , 
and  that  Eq.  (25)  was  satisfied  at  that  time.  We  can  then  push  the  ions  to  new  velocities  v 
at  and  positions  (r*,z*)  at  t'"^*,  using  the  potentials  known  at  time  T.  From  the  new  ion 
positions,  we  can  calculate  a  provisional  value  of  AN^,  which  we  shall  call  ANjjj.  AN^ 
is  slowly  varying  and  can  be  evaluated  at  the  beginning  of  the  time  step.  However,  in 
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general  Eq.  (25)  will  not  be  exactly  satisfied  at  the  end  of  the  new  time  step  (again  the  error 
will  be  of  order  Atj^),  because  the  ion  and  electron  distributions  in  phase  space  will  have 
changed  a  little  since  the  previous  time  step,  and  the  flux  condition  (12)  will  have  been  in 
error  by  order  Atj. 

This  is  indicative  that  an  increment  8  should  have  been  added  to  the  sheath 
potentials  so  as  to  preserve  the  flux  condition  throughout  the  time  step  and  the 
quasineutrality  condition  at  the  end  of  the  time  step.  In  accordance  with  the  reasoning  that 
led  to  Eq.  (22),  this  leads  to  a  correction  6Njjj  to  AN^y,  given  by 

6Nijj  =  -  (aj+aj.,)80y  +  (26) 

where  the  coefficients  aj  are  again  given  by  Eq.  (23).  Now  Eq.  (25)  can  be  written  as 

The  RHS  of  Eq.  (27)  can  be  regarded  as  known,  so  Eq.  (27)  is  an  ordinary  difference 
equation  that  determines  the  increments  to  the  sheath  potential  8^gj.  The  equation  is 
strongly  nonlinear  through  the  term  ANe?**'.  Once  the  8^y  are  known,  the  ions  are  pushed 
again  in  the  newly  determined  potentials,  and  the  electrons  with  energy  exceeding  the 
sheath  potential  are  allowed  to  escape  to  the  wall,  as  discussed  in  Sec.  4.  Thus  the 
computation  can  proceed,  with  the  ion  dynamics,  the  sheath  potentials,  and  the  transverse 
electric  fields  computed  self-consistently  at  each  time  step. 

The  reader  may  have  noticed  that  the  calculation  of  the  sheath  potential  correction 
6^sj  is  performed  fully  implicitly  in  conjunction  with  the  particle  push.  This  is  in  fact 
necessary  for  numerical  stability:  as  will  be  seen  in  the  next  section,  the  sheath  potential 
exerts  a  diffusive  effect  on  transverse  ion  current,  and  the  natural  time  step  for  this  diffusion 
is  shorter  than  we  wish  to  resolve.  Everything  else  in  the  computation  is  done  explicitly. 
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7.  Nature  of  ion  transport 


A.  Insulating  walls:  field-aligned  flow  model 

Section  6A  makes  it  clear  that  net  ion  transport  across  field  lines  is  impossible  in  an 
insulating  vessel,  since  the  plasma  must  remain  quasineutral  and  electron  transport  across 
field  lines  is  negligible.  However,  it  is  quite  possible  for  there  to  be  a  non-zero  transverse 
ion  flux  at  particular  locations  along  a  field  line.  Any  transverse  flux  at  one  location  must 
be  canceled  out  by  a  reverse  flux  elsewhere  along  the  field  line.  Thus,  in  situations  where 
there  are  strong  variations  along  a  field  line  (as  in  Fig.  2  where  the  field  lines  flare  outward 
downstream)  one  may  find  patterns  of  cross-field  ion  eddy  flow.  In  the  next  section  we 
shall  see  an  example  of  this  type  of  behavior  in  our  simulations. 

Nevertheless,  we  shall  consider  here  a  simple  solvable  one-dimensional  fluid 
approximation,  in  which  the  ion  flow  is  assumed  to  be  field-aligned.  This  model,  which 
extends  the  work  of  Godyak  and  collaborators^'  to  the  case  of  flow  along  curved  magnetic 
field  lines,  helps  to  interpret  the  two-dimensional  simulations,  and  provides  insight  into  the 
way  in  which  the  plasma  temperature,  density,  and  other  properties  are  determined  by  the 
quasineutrality  condition,  the  ambipolar  potential,  the  Bohm  boundary  condition,  and  the 
requirement  of  global  power  balance.  In  the  model,  the  electrons  on  any  single  field  line  are 
represented  as  an  isothermal  fluid.  The  ions  are  treated  as  a  cold  fluid,  with  field-aligned 
flow  velocity  Uj,  subject  to  an  electric  field  En  =  -  9^ /  9  $  and  a  constant  mean  free  path  Xj 
for  (predominantly  charge  exchange^'*)  collisions  with  stationary  neutrals.  In  addition,  we 
assume  that  electron-impact  ionization  proceeds  at  frequency  Vj,  a  function  of  T^.  In  steady 
state,  the  ion  continuity  and  momentum  conservation  equations,  and  the  Boltzmann  relation 

e[0 -^(  5=0)]  =T,ln[n/n(  5=0)1,  (28) 

can  be  reduced  to  a  single  ordinary  differential  equation  for  the  normalized  ion  flow 
velocity  u  =  u/c^. 


21 


du  I  u  I 


d  In  B 
u  — — 


V 


(29) 


ds 


X 


and  a  quadrature 

n(s) 

n(0) 


for  the  normalized  density, 

-S 


exp 


-  f  ds '  u 

n 


lul 


du  ’ 
ds ' , 


(30) 


In  Eqs.  (29)  and  (30),  v=S„,,v,/c„  and  X=Xi/$„,a„  where  $  is  the  curvilinear 

coodinate  along  the  field  line  and  <;  =0,  $  =  are  the  points  at  which  the  field  line 
intersects  the  vessel  walls. 

The  boundary  conditions  for  Eq.  (29)  arise  from  the  Bohm  condition,  but  some 
discussion  is  needed  as  to  the  proper  way  to  apply  the  Bohm  condition  within  a  model  of  1- 
D  ion  flow  along  field  lines.  If  the  field  line  intersects  the  walls  normally,  the  appropriate 
boundary  conditions  on  the  field-aligned  velocity  U;  are  given  directly  by  the  Bohm 
condition  as  stated  in  Eqs.  (16), 


Ui(0)  =  -Cs, 

But  if  the  field  line  intersects  the  walls  obliquely,  at  angles  Oq  and  at  $  =  0  and  the 
appropriate  boundary  conditions  at  the  sheath  edge  are  Eqs.  (16)  for  the  velocity  component 
normal  to  the  wall.  To  the  extent  that  flow  within  the  quasineutral  plasma  is  strictly  field- 
aligned,  this  would  then  give 

Ui(0)  =  -Cs/sin  Oq, 

^i(^max)  ~  +Cj/sin  ^niax" 

However,  Eqs.  (32)  are  clearly  not  acceptable  boundary  conditions,  since  Eq.  (29)  becomes 
singular  when  lul  =1.  The  resolution  to  this  dilemma  is  that  one-dimensional  flow  along 
field  lines  is  not  a  consistent  model,  even  within  the  quasineutral  plasma,  at  a  point  where 
the  field  lines  intersect  a  wall  obliquely.  There  is  always  a  pre-sheath  structure, beginning 


(32a) 

(32b) 


(31a) 

(31b) 


22 


about  where  Eqs.  (3 1)  are  satisfied,  wherein  the  ions  are  accelerated  toward  the  wall  and  the 
density  falls  off  by  a  factor  of  order  sin  6.  Thus  the  consistent  boundary  conditions,  within 
the  region  where  one  can  use  a  model  of  field-aligned  ion  flow,  are  indeed  Bqs.  (31). 

The  two  boundary  conditions  (31)  applied  to  the  single  first-order  ordinary 
differential  equation  (29)  are  an  indication  that  Eqs.  (29)  and  (31)  together  constitute  a 
nonlinear  eigenvalue  problem:  the  ionization  coefficient  Vj  must  take  the  value  necessary  to 
permit  a  solution.  Since  v,  is  a  rapidly  increasing  function  of  T^,  Eqs.  (29)  and  (31)  actually 
determine  the  discharge  temperature,  as  well  as  the  velocity  profile  Uj(S).  The  normalized 
density  profile  n(S)/n(0)  and  the  normalized  potential  [^(5)  -  ^(0)]  /T^  are  then  determined 
by  Eqs.  (30)  and  (28).  The  absolute  plasma  density  is  determined  by  the  power  balance 

Win  =  Wout  =  [n(0)C3Aj^(0)+n{<:^^)CsAj_(5^)  ]  [ei+e^s+-Tej  ,  (33) 

where  Wi„  is  the  microwave  power  input  to  the  flux  tube,  W<,u,  is  the  energy  loss  rate  due  to 
electrons  and  ions  reaching  the  walls,  and  A^di)  is  the  cross-sectional  area  of  the  flux  tube 
at  $.  For  Maxwellian  electrons  and  cold  ions,  the  energy  loss  per  escaped  ion  is  the  sum  of 
e,,  the  electron  inelastic  collison  energy  loss  per  ionization  (including  excitation, 
dissociation,  etc.  as  well  as  ionization);  e(P^+TJ2,  the  ion  energy  at  the  wall;  and  2Te,  the 
mean  thermal  energy  of  an  escaping  electron. 

Equations  (28)-(31)  show  that  the  normalized  profiles  Ui(<:)/Cs,  n($)/n(0)  and  [^($)- 
^(0)]/Te  depend  only  on  the  collisionality  parameter  X  and  the  shape  of  the  magnetic  field 
lines  (but  not  the  magnitude  of  the  magnetic  field).  The  temperature  T^  is  determined  by  the 
requirement  that  the  ionization  rate  be  just  sufficient  to  replace  the  ions  lost  to  the  walls. 
Since  the  eigenvalue  v  is  the  ionization  rate  Vj  is  inversely  dependent  on  the  field 

line  length  1“  “lot®  physical  terms,  ions  escape  more  rapidly  from  shorter  field  lines, 
and  thus  the  ionization  rate  (i.e.  T^)  must  be  larger  to  replace  the  lost  ions.  The  density  must 
be  lower  on  these  field  lines  so  that  a  specified  level  of  microwave  power,  divided  among 
fewer  electrons,  is  sufficient  to  raise  T^  to  the  required  value. 

In  Fig.  3,  we  show  a  solution  of  Eqs.  (28)-(31)  for  a  case  with  X  =  0.27,  $niax“3^ 

and 
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B  =  B^exp 


(34) 


—Jl—J 

(26.9  cm)'‘J  ' 

corresponding  roughly  to  one  of  the  long  interior  field  lines  of  Fig.  2,  with  the  gas  being 
argon  at  pressure  1  mTorr.  A  strong  left-right  asymmetry  is  evident,  introduced  by  the  non- 
uniform  magnetic  field.  At  z  =  0,  the  ions  are  accelerated  in  a  standard  presheath  of 
thickness  ~Xi,  but  to  the  right  the  ions  are  accelerated  by  an  electrostatic  force  arising  from 
the  mirror  effect,  and  the  plasma  density  falls  off  as  the  field  lines  diverge.  Thus  there  is 
less  need  on  the  right  for  presheath  acceleration  to  satisfy  the  Bohm  condition,  and  indeed 
the  presheath  is  less  evident. 

In  Fig.  4,  we  show  how  the  electron  temperature  on  a  particular  field  line  depends 
on  the  length  of  that  field  line,  with  B(S)  given  by  (34).  Referring  to  the  geometry 
shown  in  Fig.  2,  we  note  that  5 max  is  essentially  the  same  for  all  the  interior  field  lines  that 
terminate  on  the  end  wall,  but  decreases  steadily  as  we  go  to  the  outer  field  lines  that 
terminate  on  the  radial  wall.  Thus  the  1-D  model  predicts  that  T^  will  be  constant  on  the 
interior  field  lines,  but  will  steadily  increase  on  the  outer  field  lines.  Similarly,  the  average 
density  <  n  >  on  a  particular  field  line  is  predicted  to  be  constant  on  the  interior  lines  but 
steadily  decreasing  on  the  outer  field  lines  that  terminate  on  the  radial  wall.  The  1-D  model 
of  Eqs.  (28)  -  (31),  (33)  has  no  dependence  on  the  angles  6q  and  so  there  is  no 
discontinuity  in  plasma  properties  between  the  last  field  line  that  terminates  (nearly 
normally)  on  the  end  wall  and  the  first  that  terminates  (at  glancing  incidence)  on  the  radial 
wall.  As  we  shall  see  in  Sec.  8,  this  is  not  a  realistic  representation  of  the  2-D  flow.  In  the 
2-D  simulation,  field  lines  which  intersect  the  radial  wall  lose  ions  much  more  rapidly,  and 
thus  Tg  is  larger  and  the  average  density  < n>  is  smaller  for  field  lines  that  terminate  on  the 
radial  wall. 

We  note  that  the  rate  at  which  ions  escape  from  a  field  line  is  proportional  to  the  ratio 
of  the  density  at  the  walls  to  the  mean  density  on  the  field  line,  »}=[nj(0)+nj($max)]  ^ 

For  this  reason,  T^  and  <  nj  >  are  sensitive  to  the  density  profile  nj(  t:)/nj(0)  along  the  field 
line.  Within  our  1-D  model,  tj  is  determined  by  the  solution  of  Eqs.  (29)  -  (31),  but  in  our 
2-D  simulations  other  effects  enter  to  change  the  density  profile,  and  thereby  also  the  mean 
values  Tgj  and  <  Uj  >. 
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g  Cross-field  transport  in  grounded  conducting  vessels 


In  the  case  of  grounded  conducting  walls,  Simon^'*'^^  showed  in  1955  that  ions  can 
diffuse  across  magnetic  field  lines  at  a  rate  characterized  by  the  ion  gyroradius  and  mean 
free  path,  without  dragging  the  electrons  across  field  lines.  Quasineutrality  is  maintained  by 
electron  flow  to  the  walls  along  field  lines.  These  calculations  have  been  refined  over  the 
years,^*  and  in  the  present  context  of  weakly  magnetized  ions,  and  strong  electric  fields 
required  to  satisfy  the  Bohm  condition,  correspond  to  mobility-limited  or  inertially  limited 
transverse  ion  flow.^"  However,  the  effect  of  the  spatial  variation  of  the  sheath  potential, 
induced  by  the  ion  flow,  has  traditionally  been  omitted  from  these  models.  Our  quasineutral 
formulation  of  the  sheath-coupled  transverse  flow  lends  itself  to  an  exploration  of  this 
effect,  which  turns  out  to  be  dominant  in  the  case  of  internal  perturbations.  The  sheath 
potential  strongly  inhibits  non-uniformities  in  the  divergence  of  the  transverse  ion  flow,  as 
can  be  seen  from  the  following  considerations.  Suppose  that  momentarily  the  flow  is  such 
as  to  increase  the  number  of  ions  on  field  line  j,  more  rapidly  than  the  number  of  ions  on 
adjacent  field  lines.  To  maintain  global  quasineutrality  on  field  line  j,  it  is  necessary  to 
reduce  the  electron  flow  to  the  walls  from  field  line  j,  in  accordance  with  Eq.  (25). 
Therefore,  the  sheath  potential  on  field  line  j  increases.  But  this  also  increases  the  plasma 
potential  everywhere  along  field  line  j,  and  thus  opposes  the  net  ion  transverse  flow  onto  the 
field  line.  To  gain  analytic  insight  into  this  effect,  we  shall  consider  some  very  simplified 
model  problems  chosen  to  elucidate  particular  aspects  of  the  problem. 


Linearized  fluid  treatment 


Consider  a  plasma  contained  within  grounded  walls  at  z  =  0  and  z  =  L,  but  unbounded 
in  the  X  and  y  directions.  Let  the  magnetic  field  B,  be  uniform  and  along  the  z-direction. 
Since  we  wish  to  focus  on  the  effect  of  sheath  potentials,  which  are  normally  very  large 
compared  to  potential  variations  within  the  bulk  plasma,  we  shall  assume  the  bulk  plasma  in 
equilibrium  is  uniform  with  density  n„  plasma  potential  (equal  to  the  sheath  potential), 
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and  temperatures  and  Tj.  However,  we  recognize  the  existence  of  a  presheath  at  the  walls 
by  allowing  the  density  at  the  walls,  n^,  to  be  smaller  than  n„  by  a  factor  r?.  We  then 
consider  a  small  deviation  from  equilibrium,  with  perturbations  n(x),  ^(x)  and  u^(x)  to  the 
density,  sheath  potential  and  ion  flow  velocity,  each  of  the  form  e*^".  Assume  for 
convenience  that  the  perturbation  maintains  the  isothermal  character  of  the  plasma,  and  also 
the  ratio  r\  of  wall  to  bulk  density. 

We  begin  with  the  linearized  ion  momentum  equation,  which  takes  the  form 


n. 


9t 


ikP^ 

m. 


nii 


-ik^  +  Vi'n^u^  =  0, 


where 


Vi 


(35) 


(36) 


Vj  is  the  ion-neutral  momentum  transfer  collision  frequency,  and  the  second  term  in  (36) 
results  from  the  magnetization  of  the  ions,  i.e.  from  the  VyXB^  force.  (This  form  is 
appropriate  when  9Uy/9t«Vi.) 

We  can  relate  ^  to  n  by  linearizing  Eq.  (12).  Defining  a  normalized  distribution 
F^(vJ  by  fe(v^)  =  (no+n)Fe(v^),  the  pertubation  to  the  total  electron  flux  to  the  two  walls  is 


JOUt 

e 


2Tjn 


m. 


e<j> . 


The  perturbation  to  the  ion  flux  to  the  two  walls  is 


(37) 


j  Bohm  ^  27jnc,, 


while  the  effect  of  the  transverse  ion  flow  is  contained  in  the  last  term  of  Eq.  (12), 
Xo^dzV  •  =  ikLnoU^.  Substituting  these  quantities  in  Eq.  (12),  we  find 
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e<t)  = 


ikLme 

27?Fe(e«J 


(38) 


u 


X  • 


For  the  case  of  a  Maxwellian  distribution,  where  equilibrium  requires 


Eq. (39) reduces  to 

=-(ikL/2r7)(m,T,)‘'2u,.  (40) 

To  complete  the  calculation,  we  use  the  isothermal  relation  P;  =  nTj  and  the  perturbed  ion 
continuity  equation, 


9n 

9t 


+ 


L 


Vj-n  =  0, 


(41) 


where  v,  is  the  ionization  rate.  The  last  two  terms  of  (41)  cancel,  since  in  equilibrium 
ijn^c,  =  v,noL.  Using  (40)  and  (41)  in  (35),  we  arrive  at  an  equation  for  u^, 


92n  Ic^T.  k^LCg  9u^ 

27J  9t 


+  V. 


9t 


=  0. 


(42) 


Each  of  the  last  three  terms  of  Eq.  (42)  leads  to  a  distinct  type  of  response.  The  first 
term,  in  combination  with  the  second  term,  leads  to  sound  waves  among  the  ions.  These  are 
ordinary  sound  waves  associated  with  the  ion  pressure  and  propagating  at  the  ion  thermal 
speed,  not  the  ion  sound  speed.  The  last  term  of  Eq.  (42)  represents  ordinary  collisional 
damping.  The  third  term  of  (42),  which  arises  from  the  response  of  the  sheath  potential,  is 
the  most  interesting  term.  Taken  together  with  the  first  term,  the  third  term  leads  to 
diffusion  of  u^,  i.e.  it  is  formally  similar  to  a  bulk  viscosity  term.  The  sheath  potential  thus 
leads  to  a  dissipative  effect  which  drives  the  flow  toward  uniformity  of  the  velocity  gradient 
9u,(/9x.  In  typical  situations,  this  is  the  dominant  term.  To  demonstrate  this  formally, 
consider  the  normal  modes  of  Eq.  (42), 
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u,(x,t)  =  u,(0)e‘'^>‘"P‘, 


(43) 


where 


The  modes  are  always  damped,  and  if 


(44) 


4  (45) 

V  ItJ  k  ItJ 

they  are  purely  damped  with  no  oscillation.  Since  T^/Tj »  1,  inequality  (45)  will  normally 
be  satisfied.  For  short  wavelength  modes,  the  dissipative  effect  arising  from  the  electron 
sheath  potential  response  easily  dominates  over  the  effect  of  ion  pressure,  and  quickly 
drives  the  system  back  to  uniformity  in  transverse  velocity. 

In  a  simulation  code,  the  time  scale  for  this  diffusive  process,  273/k  Lx^,  can  be 
very  short.  In  typical  applications  of  our  code,  we  may  have  L  “  35  cm,  c^  “  3^10  cm/s, 
7j «  0.3,  and  the  shortest  wavelength  modes  have  ir/k  equal  to  the  spacing  of  the  field  lines 
used  in  the  grid,  as  small  as  0.3  cm.  This  gives  =  5x  10  '°.  Since  we  wish  to  use  time 
steps  much  longer  than  this,  it  is  necessary  to  calculate  the  sheath  potential  implicitly,  as 
discussed  in  Sec.  6. 


Nonuniform  microwave  heating 


It  is  important  to  note  that  the  dissipative  effect  discussed  in  the  preceding  paragraph 
drives  the  system  toward  uniform  transverse  flux,  not  toward  uniform  density.  This 
distinction  has  important  consequences.  Consider,  for  example,  the  "high  mode  of  an  ECR 
discharge,^’  in  which  all  of  the  injected  microwave  energy  is  absorbed.  The  power  density 
deposited  into  plasma  heating  thus  depends  only  on  the  microwave  flux,  not  on  the  plasma 
density  or  other  details  of  the  plasma  state.  The  rate  of  creation  of  new  electron-ion  pairs  by 
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electron-inipact  ionization  is  then  simply  proportional  to  the  microwave  flux.  Now  suppose 
there  is  a  small  sinusoidal  perturbation  to  the  microwave  flux,  which  leads  to  an  ionization 
rate  of  the  specified  form  njo  +  nje““.  Then  the  plasma  density  will  have  a  similar  perturbed 
form,  n^  +  ne'''*.  If  there  were  no  ion  transport  across  field  lines,  the  density  perturbation 
would  simply  be  proportional  to  the  ionization  source  perturbation. 


£  _  {46a) 

^10 

If,  on  the  other  hand,  the  ions  can  freely  diffuse  or  flow  across  field  lines,  the  non¬ 
uniformity  in  density  will  be  smoothed  out,  i.e. 


n 

no 


0. 


(46b) 


To  determine  n/iiQ  correctly,  we  use  Eqs.  (35),  (40),  and  the  linearized  continuity  equation, 


3n  . ,  TjnCg 

—  +  ikn„u,,  + - 

at  °  ^  h 


n-r 


0. 


(47) 


[Note  the  difference  in  the  last  term  on  the  LHS  of  (47),  as  compared  to  Eq.  (41)  where  the 
ionization  source  is  assumed  to  be  isothermal.]  Neglecting  the  effect  of  ion  collisions,  the 
steady  state  solution  of  Eqs.  (35),  (40),  and  (47)  is  given  by^* 


n  _  Te  (48) 

no  •  T^+2Ti- 

^10 

Comparing  with  Eq.  (46a),  we  see  that  transverse  ion  transport  does  reduce  the  density 
perturbation,  but  only  by  a  factor  T^/  (Tg+2Tj)  which  is  very  close  to  unity.  Diffusion 
toward  uniformity  in  density  is  strongly  inhibited  by  the  sheath  potential  response. 
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Ion  transport  to  vessel  walls 


We  have  seen  that  the  response  of  the  sheath  potential  strongly  suppresses  internal 
cross-field  flows  within  the  plasma.  However,  it  has  been  known  for  decades^'*  that  loss  of 
ions  to  the  walls,  through  non-ambipolar  cross-field  flow,  is  an  important  process  in 
conducting  vessels.  The  sheath  potential  response  can  influence  this  flow,  but  it  does  not 
change  the  general  conclusion.  The  effect  of  the  sheath  potential  is  merely  to  induce  a 
smooth  flow  pattern  which  has  uniform  divergence,  i.e.  reduces  the  density  at  a  comparable 
rate  on  different  field  lines. 

A  more  systematic  discussion  of  cross-field  transport,  for  discharges  within 
conducting  vessels,  will  be  provided  in  future  publications. 

8.  Simulation  of  an  ECR  discharge  in  an  insulating  vessel 


An  axisymmetric  simulation  code  QUASI-rz  has  been  developed  which  implements 
the  formalism  presented  in  Secs.  1-6.  Electrons  are  represented  as  particles  transported 
along  the  field  lines,  as  discussed  in  Sec.  3.  Each  electron  is  characterized  by  its  position  z, 
its  parallel  velocity,  and  the  value  of  its  magnetic  moment  (which  is  equivalent  to  knowing 
the  magnitude  of  its  perpendicular  velocity).  Ions  are  represented  as  particles  characterized 
by  position  in  the  r-z  plane  and  by  all  three  components  of  velocity,  as  discussed  in  Sec.  6. 
Within  the  plasma,  parallel  electric  fields  are  calculated  as  in  Sec.  3,  and  transverse  fields  as 
in  Sec.  6.  Sheaths  at  passive  surfaces  are  represented  as  thin  potential  barriers,  as  in  Sec.  4, 
while  the  Bohm  boundary  condition  is  imposed  as  in  Sec.  5.  Electron-neutral  and  ion- 
neutral  collisions  are  included,  via  a  Monte  Carlo  step  (using  the  null-collision  method^^  '^") 
which  occurs  at  the  end  of  each  particle  push  step.  Electron-electron  collisions  are  included 
in  the  Monte  Carlo  step  via  the  recently-developed  Langevin  formulation  of  multiple  small- 
angle  scattering,^'  which  is  a  great  advance  in  efficiency  over  previous  numerical 
formulations  in  terms  of  binary  collisions. 

The  code  is  ultimately  intended  to  provide  a  complete  kinetic  picture  of  an  ECR 
discharge.  We  shall  present  here  a  sample  calculation,  which  somewhat  oversimplifies  the 
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representation  of  ECR  experiments  which  are  ongoing  at  nonetheless 

illustrates  a  number  of  interesting  points.  The  geometry  is  taken  to  be  a  cylindrical 
insulating  vessel  of  length  35  cm  and  radius  7  cm,  representing  the  source  region  of  the 
ECR,  as  illustrated  in  Fig.  2.  The  magnetic  field  configuration,  as  shown  in  Fig.  2,  is  the 
actual  field  applied  in  the  experiment,  ranging  from  1  kG  at  z=0  to  171  G  at  z  =  35  cm.  The 
gas  is  argon  at  pressure  1  mTorr,  and  a  collision  set  is  used  which  includes  accurate  cross 
sections  for  electron-neutral'*^  and  ion-neutral  elastic  scattering, electron-impact  ionization 
and  excitation,'*^  and  ion-neutral  charge  exchange.^o-^*  Leonhardt,  et  al'*’  have  recently 
shown  that  interactions  between  metastables  and  other  components  are  not  important  in 
high-density  Ar  plasma  at  these  pressures,  so  these  are  not  included  in  the  collision  set. 

The  representation  of  electron  cyclotron  heating  is  simplified  in  the  present 
simulation.  We  do  not  calculate  the  propagation  of  the  microwaves  or  the  details  of  their 
interaction  with  the  electrons.  Rather,  we  heat  the  electrons  by  giving  each  electron  a 
random  kick  in  transverse  velocity  every  time  it  crosses  the  resonant  surface,  which  lies  in 
the  plane  z  =  4  cm.  The  magnitude  of  the  kick  Avj^  is  chosen  randomly  from  a  Maxwellian 
distribution  (2tr)-‘'2exp[-(Avj)2/2A2].  The  root-mean-square  value  A  is  chosen  so  that  the 
total  power  absorbed  by  all  the  electrons  on  field  line  j  is  equal  to  P^(rj,res)'^,res’  where  P^(r) 
is  the  specified  microwave  power  density,  rj  ^  is  the  value  of  r  where  field  line  j  intersects 
the  resonant  surface,  and  is  the  area  of  flux  tube  j  on  the  resonant  surface.  In  the 
present  simulation,  we  use 

P^(r)  =  Poexp(-r2/ro2),  (4^) 

with  Po  =  1.8  W/cm^  and  r^  =  7  cm.  The  total  absorbed  microwave  power,  integrated  from 
r = 0  to  the  wall  at  r =7  cm,  is  then  350  W. 

Given  the  simplication  of  the  geometry  and  of  the  ECR  heating  process,  one  cannot 
necessarily  expect  the  simulation  to  provide  a  quantitatively  acccurate  picture  of  any 
particular  experiment.  However,  the  model  includes  all  of  the  important  effects  that 
determine  the  plasma  properties,  including  spatial  variation  of  parameters,  flow  patterns, 
velocity  distribution  functions,  and  ionization  fractions. 
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Simulation  results 


We  chose  (arbitrarily)  to  initiate  the  simulation  with  uniform  temperature  =  6.6  eV 
and  plasma  density  n(r,z)  proportional  to  B(z,r),  with  n  =  1 10*^  at  z  =  0.  After  a  time  of 
order  300  ns  the  simulation  evolves  to  a  steady  state,  which  appears  to  be  insensitive  to  the 
initial  conditions.  We  shall  show  graphics  illustrating  the  steady  state  plasma  at  time  500 
ns. 

Figure  5  shows  the  plasma  density  n(r,z).  The  peak  density  is  6.3x10*^  cm■^  in 
reasonable  agreement  with  experimental  results.  Curiously,  the  peak  is  seen  to  be  off-axis, 
with  a  20%  density  dip  on  axis,  even  though  the  microwave  power  is  gently  peaked  on-axis. 
As  we  shall  see,  this  is  a  consequence  of  two-dimensional  flows  which  differ  from  the  field- 
aligned  flow  model  of  Sec.  7A. 

Figure  6  shows  surface  plots  of  the  ion  fluid  velocity  components  iiz(r,z)  and  UiXr,z). 
We  note  that  ions  flow  to  all  walls  at  flow  velocity  U;^  “  c^,  in  accordance  with  the  Bohm 
condition.  At  every  wall,  a  presheath  is  evident  in  Figs.  5  and  6,  wherein  the  ions  are 
accelerated  to  c^,  but  the  presheath  is  minimal  on  the  downstream  end  wall  where  the 
acceleration  largely  occurs  as  a  consequence  of  magnetic  field  expansion.  Figure  7  is  a 
vector  plot  of  the  ion  flux,  which  more  clearly  illustrates  the  surprisingly  complex  nature  of 
the  flow.  The  important  feature  is  that  the  outer  field  lines  of  Fig.  2,  which  intersect  the 
radial  wall  at  an  acute  angle  of  5  *  to  19* ,  lose  ions  to  the  wall  over  a  relatively  large  area  at 
perpendicular  velocity  Uj"*" = Cj.  This  loss  of  ions  is  not  resupplied  through  flow  along  the 
field  line,  as  in  the  model  of  Sec.  7A,  but  mainly  through  cross-field  flow  from  adjacent 
field  lines.  However,  the  quasineutrality  condition  specifies  that  no  field  line  can  have  a  net 
gain  or  loss  of  ions  through  cross-field  flow  within  the  plasma.  Thus,  there  must  be  an 
inward  return  flow  upstream.  Comparing  Fig.  7  with  Fig.  2,  it  can  be  seen  that  the  result  is 
a  flow  pattern  which  is  fairly  close  to  field  aligned,  and  yet  shows  a  significant  cross-field 
eddy  flow  structure. 

The  nature  of  the  flow  has  a  significant  impact  on  the  density  and  temperature  profile 
within  the  plasma.  The  cross-field  flow  permits  the  rate  of  ion  loss  from  those  field  lines 
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that  terminate  obliquely  on  the  radial  wall  to  be  greater  than  the  loss  rate  from  the  interior 
field  lines  that  terminate  nearly  normally  on  the  end  wall.  (Recall  that  in  the  1-D  model  of 
Sec.  7A,  the  loss  rate  depends  only  on  the  length  of  the  field  line,  but  not  on  the  angle  at 
which  the  field  line  intersects  the  wall.)  Thus  there  is  an  abrupt  decrease  in  the  plasma 
density,  as  we  move  from  the  last  field  line  that  terminates  on  the  end  wall  (j  —  9)  to  the  first 
field  line  that  terminates  on  the  radial  wall.  This  is  evident  in  Fig.  5.  In  Fig.  8  we  plot  the 
average  value  of  the  density  <nj>  on  field  line  j,  as  a  function  of  j.  This  figure  shows  even 
more  clearly  the  sharp  change  in  the  plasma  around  field  line  j  =  9.  As  j  increases  further, 
i.e.  as  we  consider  the  outer  field  lines  that  terminate  at  decreasing  values  of  z,  ions  are  lost 

still  faster  and  there  are  further  decreases  in  <  nj  >. 

The  electron  temperature  T^,  which  is  found  to  be  quite  close  to  isotropic,  is  shown  as 
a  surface  plot  in  Fig.  9.  T^  is  close  to  flat,  as  expected,  on  the  inner  field  lines  j  ^  9.  On 
field  lines  j  >  9,  which  terminate  on  the  radial  wall,  T^  increases  as  the  field  line  becomes 
shorter.  This  is  a  consequence  of  the  more  rapid  loss  of  ions  on  these  short  field  lines,  as 
discussed  in  Sec.  7A.  However,  the  electron  temperature,  which  we  define  as  the  mean 
<3meV^/2>,  is  found  to  be  significantly  larger  than  the  values  obtained  from  the  fluid  model 
of  Sec.  7A.  The  primary  reason  for  this  is  that  the  electron  energy  distribution  function  falls 
off  from  Maxwellian  in  the  high  energy  regime,  which  controls  electron-impact  ionization. 
Thus  Tg  must  be  somewhat  larger  to  give  the  appropriate  ionization  rate.  Distribution 
functions  will  be  discussed  in  more  detail  in  forthcoming  work.  This  temperature  and 

density  pattern  is  seen  in  the  NRL  experiments."*^ 

The  cross-field  ion  flow  also  leads  to  a  rather  subtle  effect  on  the  density  of  the 
interior  field  lines.  At  large  z  (downstream),  the  flow  is  divergent  and  depletes  the  ions,  but 
in  the  region  of  highest  density  near  z  =  10  cm,  the  flow  is  convergent  and  brings  in 
additional  ions.  Thus  this  flow  tends  to  increase  the  peakedness  of  the  density  profile  nj(z) 
along  a  field  line,  i.e.  to  increase  the  maximum  of  Uj  but  reduce  nj  at  the  ends  of  the  field 
line.  Figure  7  shows  that  the  cross-field  eddy  flow  is  strongest  on  the  field  lines  (j=7,8,9) 
that  terminate  at  the  outside  of  the  end  wall.  On  field  lines  that  are  near  the  axis,  the  flow  is 
very  nearly  field-aligned.  Thus  the  density  profile  nj(z)  is  more  peaked  on  field  lines 
j=7,8,9  than  on  the  inner  field  lines  near  the  axis.  However,  the  loss  rate  of  ions  to  the  walls 
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is  proportional  to  the  density  at  the  end  of  the  field  line,  and  thus  is  slower  on  field  lines 
j=7,8,9.  Applying  the  arguments  of  Sec.  7A  in  this  more  general  multidimensional  context, 
reducing  the  particle  loss  rate  increases  the  density  all  along  the  field  line.  This  is  the 
explanation  of  the  off-axis  density  peak  which  can  be  seen  in  Fig.  5  at  z  =  8,  r  —  2.8,  and 
which  is  also  evident  in  Fig.  8. 

The  potential  profile  0(r,z),  shown  in  Fig.  10,  is  more  complex  than  might  have  been 
anticipated.  It  is  not  peaked  on  axis,  and  it  is  not  simply  indicative  of  the  electron  pressure 
profile,  as  in  an  unmagnetized  plasma.  Rather,  <f>{v,z)  has  a  saddle  at  on  axis  at  z  =  9,  and  for 
z<  15  cm  ^  is  a  monotonically  increasing  function  of  r,  except  in  the  presheath  at  the  radial 
wall.  Downstream,  0  is  a  monotonically  decreasing  function  of  r.  This  structure  is  just 
what  is  needed  to  drive  the  radially  outward  ion  flow  downstream,  and  the  inward  return 
flow  upstream.  This  type  of  potential  structure  has  been  observed  in  the  NRL 
experiments.'^ 


9.  Concluding  Remarks 

We  have  presented  a  variety  of  techniques  for  modeling  the  quasineutral  region  of  a 
plasma  discharge.  These  include  methods  for:  (i)  determining  the  electric  field  parallel  to 
B,  and  the  associated  electron  transport;  (ii)  determining  the  potential  variations  transverse 
to  B,  and  the  associated  ion  transport;  (iii)  determining  the  spatially  and  temporally 
dependent  sheath  potential;  (iv)  enforcing  the  Bohm  flow  condition.  All  of  these 
techniques  are  time-dependent  and  fully  kinetic,  and  based  on  PIC  modeling  of  both  the 
electrons  and  ions.  In  all  cases,  the  fields  and  potentials  are  determined  directly  from  the 
requirement  of  quasineutrality.  Poisson’s  equation  is  not  used,  the  Debye  shielding  length 
is  effectively  set  to  zero,  and  electron  plasma  oscillations  are  not  present  in  the  modeling. 
As  a  result,  the  models  focus  on  the  macroscopic  processes  that  are  aetually  of  interest,  and 
permit  the  use  of  large  time  steps  and  spatial  gridding  that  enhance  computational 
efficiency. 
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The  modeling  techniques  have  been  discussed  in  the  context  of  magnetized  ECR 
discharges  for  processing  applications,  but  we  believe  that  they  lend  themselves  to  a  wide 
range  of  plasma  conditions,  including  unmagnetized  plasmas,  unbounded  plasmas,  and 
high-temperature  collisionless  plasmas  of  interest  in  fusion  and  space  physics. 

The  power  of  the  techniques  has  been  demonstrated  by  applying  them  analytically  to 
some  simple  fluid  situations,  and  by  incorporating  them  in  a  two-dimensional 
(axisymmetric)  PIC  simulation  model  for  an  ECR  discharge.  Analytic  study  of  discharge 
contained  within  a  conducting  vessel  revealed  significant  modifications  to  the  classic  Simon 
diffusion  across  a  magnetic  field.^^-^®  Simulation  of  an  ECR  discharge  contained  within 
insulating  walls  revealed  unanticipated  and  important  structural  features,  which  are 

dependent  on  2-D  ion  flows. 

In  this  paper  we  have  not  discussed  the  modeling  of  the  microwave-plasma 
interaction.  Computationally  efficient  approaches  to  this  problem  will  be  discussed  m 
subsequent  publications.  Simulation  of  discharges  within  conducting  vessels  are  also  in 
progress  and  will  be  reported  in  future  publications. 
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Table  1.  TYPICAL  PARAMETERS  FOR  AN  ECR  REACTOR 


Microwave  frequency 
Pressure 
Magnetic  field 
Plasma  density 
Electron  temperature 

Nominal  Time  Scales 


Electron  plasma  oscillation  10  ‘  ‘  sec 

Electron  gyration  6x10'*’  sec 

Electron  transit  through  ECR  zone  10  ®  sec 
Electron  collision  frequency  10  ’  sec 

Ion  plasma  period  10  ’  sec 

Ion  gyration  10'^  sec 

Ion  collision  frequency  10'^  sec 

Ion  lifetime  lO  '*  sec 

Neutral  residence  time  10  ’  sec 


f„  ~  2.45  X  10’  Hz 
P  ~  0.5  - 10  mTorr 
Bq  ~  20  -  1000  Gauss 
Ug  ~  Uj  ~  10' ‘  -  lO'Vcm’ 
Tg  ~  1-10  eV 


Spatial  Scales 


Debye  length 

<  0.01  cm 

Electron  gyroradius 

0.005  to  0.2  cm 

Sheath  thickness 

0.1  cm 

Electron  mfp 

few  to  10’ s  of  cm 

Microwave  wavelength 

<  12  cm 

Ion  gyroradius 

up  to  few  cm 

Ion  mfp 

mm’s  to  cm’s 

Device  size 

10  - 100  cm 
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Fig.  2  —  In  the  simulation  shown  in  this  paper,  the  geometry  of  the  ECR  source  is 
simplified  to  a  cylinder.  The  simulation  grid,  shown  in  this  figure,  is  formed  by  a 
chosen  set  of  magnetic  field  lines  (the  actual  measured  field  lines  in  an  experimental 
configuration),  and  by  an  equally  spaced  z-grid. 
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Fig.  3  —  Plasma  state  obtained  by  solution  of  the  1-D  field-aligned  flow  model,  Eqs.  (28)  -  (31), 
for  conditions  approximating  the  long  field  lines  of  Fig.  2.  Shown  are  the  normalized  density 
N(s)  s  n(s)/n(0)  (dashed  curve),  the  normalized  ion  flow  velocity  u(s)  s  Ui(s)/Cs  (solid  curve),  and 
the  normalized  plasma  potential  <I*(s)  =  [(j)(s)-(|)(0)]/Tg  (dot-dashed  curve).  The  eigenvalue  is  v  =  1.56, 
corresponding  to  Tg  =  3.4eV. 
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(b) 

Fig.  6  —  Surface  plot  of  the  ion  fluid  velocity  components  (a)  U;^  and 
(b)  Ujf  at  time  500  ps,  from  the  simulation  of  Fig.  5. 
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7  —  Vector  plot  of  the  ion  flux  at  time  500  psec,  from  the  simulation  of  Fig.  5 
Note  that  the  aspect  ratio  of  the  plot  is  different  from  that  of  the  simulation. 
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